% Function that calculates the frequencies of correct identification.
% Arguments:
%   M       Number of simulations
%   T0      Minimum series length.
%   T1      Maximum series length.
%   e1fun   function that controls the rate of reduction of the parameter
%           of the approximate algorithm 2sn version (the one that separates the roots)
%   e2fun   function that controls the rate of reduction of the parameter
%           of the approximate algorithm 2sn version (the one in the Smith form estimation)
%   n       dimension of the series
%   cases   DGP names
%   rts     unit roots of the DGPs: first row 1, second row -1.
%   ord     order of the DGPs
%   joh     Do Johanson test or not
%
% Output:
%   pcmat   Frequency of right identification.
%   pcvmat  Frequency of right identification of the unit roots (not structure).
%   pcomat  Frequency of right identification of the order of the model.
%   pcjmat  Frequency of right identification by Johansen.
%   Tgrid   Grid of series lengths.

function [pcmat, pcvmat, pcomat, pcjmat, Tgrid]=sequence_simul_nodiv2(M,T0,T1,e1fun,e2fun,n,s,cases,rts,ord,joh)

nc=length(cases);

hT=25;
nT=(T1-T0)/hT+1;

pcmat=zeros(nc,nT);
pcvmat=zeros(nc,nT);
pcomat=zeros(nc,nT);
pcjmat=nan(nc,nT,2);

Tgrid=linspace(T0,T1,nT);

alpha=[0.01 0.05];

for t=1:length(Tgrid)
    
    t

    for i=1:nc
        
        if i==2,
            1;
        end
                
        T=Tgrid(t);

        e_1=feval(e1fun,T);

        e_2=feval(e2fun,T);
        
        pc=0;
        pcv=0;
        pco=0;
        pcj=[0 0]';

        for j=1:M

            flag=1;

            while flag

                [A, r, p, x]=feval(cases{i},T); % get simulated data
                
                nrv=sum(1./abs(r)<1+e_2);       % number of unit roots
                if nrv==sum(sum(rts{i})),
                    pcv=pcv+1;
                end
                
                if p==ord(i),                   % Correct order
                    pco=pco+1;
                end

                [U, D0, V]=smith_polymat_nodiv(A,e_1); % First Smith form estim.
                                   
                rv=1./r;
                
                urv=rv(abs(rv)<1+e_2);                  % Get unit roots.
                
                if ~isempty(urv),                        
                    
                    urvp=maproots(urv,2);               % Project roots.
                    
                    D1=maproots_mat(D0,urvp);           % Assign roots in matrix.
                    
                    [U, D, V]=smith_polymat(D1,0.001);  % Final Smith form estim.
                    
                    % Unit roots count.
                    
                    for k=1:n
                        for l=0:s-1
                            nr(l+1,k)=sum(abs(roots(D{k,k})-exp(2*pi*1i*l/s))<1.001);
                        end
                    end
                    
                else
                   
                    nr=zeros(s,n);
                    
                end
                                
                if nr==rts{i},  % Right unit roots structura identification.
                    pc=pc+1;
                else
                    1;
                end
                
                if joh
                    
                    %%%%%% Johansen test %%%%%%
                                                           
                    [h,pValue,stat,cValue,mles] = jcitest(x','display','off');
                    
                    % We choose the cointegration rank using the rejection
                    % vector h. (h=1 -> reject, h=0 -> do not reject).
                    % Cointegration rank = maximum rank for which the test
                    % do not reject.
                                                            
                    pValueVec=nan(1,n);
                    for k=1:n
                        pValueVec(k)=pValue{1,k};
                    end
                    
                    rechaza=(ones(2,1)*pValueVec<alpha'*ones(1,n));
                    rest=sum(rechaza>0,2);
                    
                    pcj=pcj+(rest==n-sum(sum(rts{i})));
                    
                end
                
                flag=0;

            end

        end

        pc=pc/M;
        pcv=pcv/M;
        pco=pco/M;
        pcj=pcj/M;

        pcmat(i,t)=pc;
        pcvmat(i,t)=pcv;
        pcomat(i,t)=pco;
        pcjmat(i,t,:)=pcj;

    end
    

end
